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Abstract 

We consider fractional relaxation and fractional oscillation equations 
involving Erdelyi-Kober integrals. In terms of Riemann-Liouville inte¬ 
grals, the equations we analyze can be understood as equations with time- 
varying coefficients. Replacing Riemann-Liouville integrals with Erdelyi- 
Kober-type integrals in certain fractional oscillation models, we obtain some 
more general integro-differential equations. The corresponding Cauchy- 
type problems can be solved numerically, and, in some cases analytically, 
in terms of Saigo-Kilbas Mittag-Leffler functions. The numerical results 
are obtained by a treatment similar to that developed by K. Diethelm 
and N.J. Ford to solve the Bagley-Torvik equation. Novel results about 
the numerical approach to the fractional damped oscillator equation with 
time-varying coefficients are also presented. 
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1. Introduction 

In this paper, we consider fractional relaxation and fractional oscilla¬ 
tion models with time-varying coefficients, involving Erdelyi-Kober-type 
integrals. In the recent years, a number of papers have been devoted to the 
analysis of fractional relaxation and oscillations as well as to their applica¬ 
tions, see, e.g., mm and references therein. For a probabilistic approach 
to fractional relaxation equations, we refer to [2]. Non-Debye relaxation 
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phenomena in dielectrics are often modeled by means of fractional differ¬ 
ential equations that generalize the classical relaxation equation. In this 
framework, Mittag-Lefffer-type functions play a relevant role to describe 
anomalous relaxation, including special cases such as the Cole-Cole [3|, the 
Davidson-Cole [6], and the Havriliak-Negami mm models (see also m 
for a short review about analytical representations of relaxation functions 
in non-Debye processes). In recent theoretical investigations [3J [IT], the 
authors studied fractional relaxation models with time-varying coefficients 
resorting to different approaches. In the first part of this paper, we consider 
a further approach in this direction, studying a new fractional model in¬ 
volving Erdelyi-Kober-type integrals. In this framework, we provide both 
analytical and numerical results. 

In the second part of the paper, we consider a fractional generalization 
of the damped oscillator equation with time-dependent elasticity involving 
Erdelyi-Kober-type integrals. Fractional generalizations of the damped 
oscillator equation have been subject of recent research, especially in the 
framework of fractional mechanics; we refer, in particular, to the recent 
books [2T) 2S] and references therein. Moreover, we observe that a connec¬ 
tion between fractional diffusion equations involving Erdelyi-Kober-type 
operators and the density law of the generalized grey Brownian motion has 
been recently studied in [[28]. This fact stresses the role of Erdelyi-Kober 
integral operators in the future developements of fractional calculus. 

As far as we know, only few works concerning analysis as well as numer¬ 
ical treatment of fractional damped oscillators with time-dependent coeffi¬ 
cients can be found in the literature. We will show that some generalizations 
of the classical equations of mechanics by means of Erdelyi-Kober-type in¬ 
tegrals imply the analysis of integro-differential equations with time-varying 
coefficients involving Caputo-type fractional derivatives. This topic seems 
to be promising and interesting both, for numerical and analytical stud¬ 
ies of fractional differential equations with variable coefficients, and their 
applications in mechanics. 


2. Fractional relaxation models involving Erdelyi-Kober-type 

integrals 

A reasonable generalization of the classical relaxation model can be ob¬ 
tained replacing the ordinary time derivative with the the Caputo fractional 
derivative, defined by (see [23J) 


D?m 


i 


T(m — a) 


rt j m 


0 


( 1 . 1 ) 
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where a > 0, m = [a], and / e AC m [a,b\. Therefore, the so-called 
fractional relaxation equation is simply given by 

D]~ a u{t ) = -Xu(t), A > 0, t> 0, a G [0,1). (1.2) 

In this case, the fractional model describes an anomalous relaxation pro¬ 
cess with a Mittag-Leffler decay, that is with an asymptotically power law 
behavior. 

Since for absolutely continuous functions, we have 

Jl~ aD \~ a u(t) = u(t) - u( 0), (1.3) 


where 




m j o 


{t-s) 13 1 /(s) ds, 


(1.4) 


is the Riemann-Liouville fractional integral of order (3 > 0 (see, [19], e.g.), 
we can write (1.2), in the equivalent integral form, 


u(t) — u( 0) = —A j\ a u(t). 


(1.5) 


In this section we consider a generalization of the fractional relaxation 
equation in the integral form (1.5), where the Riemann-Liouville integral 


is replaced by the Erdelyi-Kober-type fractional integral (see for example 
[271 EH [32!) 

iTf = r( v / {t m - s m r~ l s m V{s) d(s m ), rj > 0, m > 0. (1.6) 
r(a) Jo 

By this, we will be able to obtain some general results concerning 
anomalous relaxation with time-dependent coefficients and hereditary ef¬ 
fects, including both, classical (ordinary) and fractional relaxation models 
as special cases. 

More general Erdelyi-Kober operators, with /3 > 0 instead of m in 
(1.6), have been studied, and the corresponding “Erdelyi-Kober fractional 
were introduced in 


m 


derivatives” were introduced in [20j, and more recently analyzed, e.g 
[ 23] . Applications of such operators in the framework of fractional mechan¬ 
ics should be the subject of further research. 

For simplicity, in (1.6) we will assume m = 1, being a G (0,1). In 
fact, in this section we consider purely relaxation phenomena. Moreover, 
we assume that the relaxation coefficient there follows a power-law behav¬ 
ior, i.e., A(t) = At 1-2 ". We are interested to understand the interplay 
between dynamical coefficients and Erdelyi-Kober-type integrals in gen¬ 
eralized mechanical models. Therefore, our first model is based on the 
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following generalization of equation (1.5) 


A 


u (t,) -U 0 = — a u(t) 


1 2 ' 

Xt-a-v rt 

r(l - at) Jo 


(t — s ) a s v u(s ) ds 


= -A 


(1.7) 


We observe that, for rj = 0 and a = 0 we recover the classical relaxation 
model. Hence, the effect of introducing Erdelyi-Kober-type integrals in 
the fractional relaxation equation is essentially that of providing power law 
time-varying coefficients. We will show that this integral equation can be 
treated in a similar way as that considered in [3j. 

Equation (1.7) can be handled in two equivalent ways. The first one is 


based on the solution of the following system of coupled integro-differential 
equations 

(u(t) - u 0 = -A t~ v ~ a D^g(t) 

I D t g = fu^t), 


( 1 . 8 ) 


which is clearly equivalent to (1.7) (here D t g = D} := dg/dt and g(0) = 0). 
This approach is interesting in view of a numerical integration of equation 
(1.7), and is suggested by a method developed by Diethelm and Ford in 


[51 [9] to solve numerically the well-known Bagley-Torvik equation (T| . 

The second approach, based on simple analytical manipulations, allows 
to reduce equation (1.7) to that considered recently in [3]. Indeed, applying 


the first order derivative to both sides of (1.7), we obtain 


d 


dt 

We will study the case 


- [u(f) + At"*-** Jt a (t v u)\ = 0. 


t TI+a u(t) + AJ t 1_a (f r? u) = 0. 


(1.9) 


( 1 . 10 ) 


Let us recall that 




( 1 . 11 ) 

and define r(t) := t r?+ "u(t). Then, taking the Caputo derivative of order 
1 — a of both sides of (1.10), we obtain the equation 

( 1 . 12 ) 


D\~ a r(t ) = -A t~ Q r{t). 


Such equation has recently been studied by Capelas de Oliveira et al. in 
[3], and its explicit solution can be written in terms of the so-called Saigo- 
Kilbas Mittag-Leffler function, for a G (0,1/2) (see the recent monograph 
m for more details). In particular, here we recall the following useful 
result (see, e.g., [T9], pp. 232-233). 
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Lemma 1.1. The function 


E 


a, l+± 


-A t a+fi 


oo k— 1 

1 + ^{-\) k t k{ - a+fi) n 

k= 1 j =0 


r ( Q (j+j| + |) + 1 ) 

r (« (i + ii + l + 1 )+ 1 )’ 

(1.13) 


solves the following fractional Cauchy problem: 

(Dfy ( t ) = —A t@y (t), t > 0, a E (0,1] , —a < fi < 1 — a, 

i 2/(0) = 1. 


(1.14) 


We recall that, according to Kilbas et al. |19] (Remark 4.8, p. 233), 
uniqueness of the solution to the Cauchy problem in (1.14) has been estab¬ 
lished only for /3 > 0. 


In view of Lemma 1.1, a solution to equation (1.12), for a E (0,1/2), is 
given by 


r(t) = E x _ a i_ « _ « (—At 1 2a ) , 
and going back to the original function, we have that 

E x _ al _^ _au (-At 1 ’ 20 ) 

1 ’ 1 —o: ’ 1-a x t 


u(t ) = 


fOt+l] 


We observe that, for a = r) = 0, equation (1.16) becomes 


u(t) = Eg i )0 (-At) = e 


—At 


(1.15) 


(1.16) 


(1.17) 


which coincides with the solution of the classical relaxation equation with 
initial condition u( 0) = 1. 

The study of the conditions for the complete monotonicity of a solution 


to the fractional Cauchy problem in (1.14) is a central topic for the physical 


sense of the anomalous relaxation models, since it ensures that in isolated 
systems the energy decays monotonically (see (T5j for a complete discussion 
about the physical realizability of the viscoelastic models). 

In [3], the authors have discussed the complete monotonicity of the 
so-called Saigo-Kilbas Mittag-Leffler function (1.13). In particular, they 


proved that the solution (1.13) is completely monotonic provided that a E 


(0,1] and /3 E (—a, 1 — a]. In our case, these conditions correspond to the 
assumption that a E (0,1/2), that ensures the complete monotonicity of 
the solution. We refer to [25] and references therein for a recent discussion 
concerning the relevance of completely monotonic functions in dielectrics. 
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1.1. Numerical results 


In this section, we present some numerical results concerning the so¬ 
lution of the fractional Cauchy problem in (1.14). In Fig. [2j we show the 


numerical solutions to such problem, for different values of a 6 (0,1], and 
varying the parameter f3 , with —a < f3 < 1 — a 

The decay of the response function, for a fixed a, is faster than expo¬ 
nential for short times, and power-law-like asymptotically, as pointed out 
also in [3J. The decay is clearly faster for lower values of a. 

In FigjTJ we observe that, for a fixed value of a (a = 0.9) for short times 
the solution to equation (1.13) decays faster with respect to the classical 


exponential relaxation decay. 

On the other hand, for a fixed value of a, the role played by the pa¬ 
rameter /3 is clear: increasing (3 the solution decays faster (see Figj2]). In 
order to understand the physical meaning of the fractional relaxation model 
described by the Cauchy problem in ( |1.14 ), we observe that, for a = 1, it 
reduces to 

\D t y (t) = -XtPy(t ), t > 0, P > -1 

(y (o) = l, (L18) 

that is to a relaxation equation with the time-dependent coefficient A (t) : = 
A t@, whose solution is given by 

( tP +1 

y(t ) = exp -A ■ 


P + 1 

From this special case, we can understand the role of the parameter (3 in 
the solution, as confirmed by the numerical results. 

In Fig.s[l]and[2j we plotted the numerical solution of (1.14) for several 
values of the parameters. The numerical method we used is an adaptive 
improvement of the predictor-corrector method earlier introduced in |7j, 
and developed in (5]. Starting from an arbitrary discretization step, h, we 
chose locally a step size inversely proportional to the size of the (classical) 
derivative of the solution being computed, so that, at the z-th step, the 
time step 

h . ■= _ ^ _ 

*' \u(ti)-u(ti- or 

i = 1,2,..., N, ti = hj, to = 0, is used. Here N is the number of the 
integration nodes, and 


Ci = 


u'{ti-2 ) - 3 ) 


U' a (ti- 1) - v!(ti-2Y 
where u’ denotes the (classical) first derivative of u ’, see [5j. 
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Figure 1. Numerical solution of (1.14) for a = 1 and a = 
0.9, A = l,/3 = 0.5, and initial discretization step h = 0.001, 
computed for t £ [0, 5]. 


In Fig. [3j we show the convergence rate for the numerical solution of 


(1.14), plotting log \u'/u\ vs j3. The graph exhibits a linear negative relation 
between the values of /3 and the rate of convergence defined as before. 


2. Damped fractional oscillator involving Erdelyi-Kober-type 

integrals 


In this section, we consider a fractional generalization of the damped os¬ 
cillator equation, involving Erdelyi-Kober-type integrals. There are many 
investigations nowadays about fractional damped oscillators involving Ca- 
puto derivatives, in the framework of the fractional mechanics (see for ex¬ 
ample [mum ini] and references therein). Here we consider a special form 
of the fractional damped equation with time-varying coefficients, contain¬ 
ing Erdelyi-Kober-type integrals, and we assume the elastic term to be 
time dependent. Hereafter, Df will denote the first order ordinary time de¬ 
rivative. A Caputo-type fractional generalization of the damped oscillator 
equation is thus given by 


D^Dtx(t) + A Dtx{t) + kx(t) = 0, (2.1) 


where A > 0 represents the viscosity coefficient and k > 0 is the elasticity 
constant. Since c Dfx(t) = J f 1- “ D t x(t), ( |2.1[ ) can also be written as 

Jl~ a (D^x) (t) + A D t x(t) + kx(t) = 0. (2.2) 


Our purpose is to consider a further, more general model of damped oscil¬ 
lator with memory effects and time-varying coefficients, replacing in (2.2) 
the Riemann-Liouville integral with the Erdelyi-Kober-type integral 

(t) = J^(t- s)" q s ?? r(s) ds. (2.3) 

Moreover, we will assume that the elasticity term depends on time, and, 
in particular, k(t) = ko/t, where we set fco := A(r/ + 1 — a) > 0. We 
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Figure 2. Numerical solution of (1.14) for a = 0.9, A = 1, 
initial discretization step h = 0.001 and several values of /3, 
computed for t e [0, 5]. 
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Figure 3. Rate of convergence, represented by \og\v!/u\. 


for the numerical solution of (1.14) for a = 0.9, A = l,/3 £ 
[—0.6,—0.1], initial discretization step h = 0.1, computed 
fro t G [0, 500], in a log — log scale. 


thus obtain the following fractional damped oscillator equation with time- 
varying coefficients, 

[DfPDt + Xt v+1 ~ a D t + kot v ~ a ) u{t) = 0. (2.4) 

Note that, for 77 = 0 and a = 1, (2.4) reduces to the classical damped oscil¬ 
lator equation with time-dependent elastic term. At best of our knowledge 
no papers can be found in the literature concerning fractional damped os¬ 


cillator equations with time-varying coefficients. Looking at (2.4), we can 


see that, replacing the Riemann-Liouville integral with the Erdelyi-Kober 
integral, power law time-variable viscosity and elasticity coefficients appear. 

The physical meaning of having such a combination (in k$) of the orig¬ 
inal constants, a and rj, is that elasticity, damping, and viscosity effects 
turn out to be coupled, that is, we may expect that damping ratio and 
(local) oscillation frequency (if any) are linked. However, we will see from 
the numerical simulations that oscillatory effects strongly depend on the 
order a of the fractional derivatives involved in the model equation. 

We now show how this choice is useful from the mathematical point of 


view. Equation (2.4) can be rewritten, decoupling it into the two auxiliary 
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coupled equations 


Dfg(t) = -A t^ +1 ~ a u(t), 
t v D t u(t) = D t g(t). 


(2.5) 


This can be shown just by taking the fractional derivative of order a of 
both sides of the second equation of the system (2.5), and replacing there 
Dfg as given by the first equation, recalling that = A(?/ + 1 — a). We 
obtain 


(2.7) 


D?D t g(t) = D t D?g(t ) = —A D t +1 ~ a u(t )) = D? (?D t u(t )), (2.6) 

where we exploited the fact that, being g{ 0 ) = 0 , ordinary and fractional 
derivatives commute. Moreover, in this case the Riemann-Liouville and 
the Caputo fractional derivatives of g(t) coincide (see, e.g., [29, pp. 73-74], 
This seems to be an interesting result, since equation (2.4) is a fractional 
differential equation with time-varying coefficients, classically related to 
fractionally damped oscillation. 

The previous equation can be written as a “hybrid” system, in the sense 
that it involves both, classical and fractional derivatives. Here, the range 
of the parameters is0<a<l,A>0, and rj > 0. Writing v(t) := D t u(t), 
we obtain 

D t u(t) = v(t) 

Df ( t 71 v) (t) = —A t s v — gt^u, 
where we set, for short, 5 := g + 1 — a, g := A (rj + 1 — a), and 7 := g — a, 
and assume as initial conditions rt(0) = u(0) = 1. Setting w(t) := 
we have 

D t u(t ) = t~ v w(t) 

Dfw(t) = —A t s ~ v w(t) — gt^u(t), 

This approach, adopted to solve numerically \2A \ , is suggested again by 
the method of Diethelm and Ford in [ 8 ] , extended according to the adaptive 
pattern developed in 0 ), to solve the celebrated Bagley-Torvik equation 

m, 

AD^y + B c D 3 t /2 y + Cy = f(t ), (2.9) 

where A ^ 0, B, and C are real constants. We recall that in [ 8 ] equation 
(2.9) was rewritten as 


( 2 . 8 ) 


I) 

I) 

D 


4/2 

t y 1 = 2/2 

1/2 

t 2/2 = yz 

\' 2 r 
R/2 _ 4-1 


( 2 . 10 ) 


D^y 4 = A-'i-Cy!- By 4 + 


with the initial conditions yi(0) = yo, 2 / 2 (0) = 0, 73 = y' Q , 2 / 4 ( 0 ) = 0. 
However, using Griiwald-Letnikov discretization, we face the problem of 
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having unevenly spaced nodes, which problem can be solved resorting to 
some interpolation. The latter can be merely linear, since the additional 
error introduced is of the second order, while the order of the method of 0 is 
of order 3/2. Quadratic interpolation was also tried, which introduces third- 
order (hence negligible) errors, at the price of some little complications. 


3. Numerical results 


( 2 . 8 ) 


In this section, we present some results for the numerical solution of 
We first discuss the role played by the parameters a, A, and ?y enter¬ 


ing equation (2.8). The fractional order of the derivative has an intrinsic 
damping effect on oscillatory phenomena, as discussed for example in [33]. 


In order to understand the role of the parameters A and 77 , we observe that 


for a = 1 , ( 2 . 8 ) becomes 


Dl + 


(? + ■') 




u(t ) = 0 , 


(3.1) 


where we recall that we assumed A > 0 and i] > 0. This is a damped 
oscillator with time-dependent viscosity and elasticity terms. We remark 
that in our fractional model, time-variable viscosity, damping, and elas¬ 
ticity terms are coupled. This implies that damping and (instantaneous) 
oscillation frequency are linked. When 77 = 0 and a = 1, (2.8) reduces to 
the classical relaxation equation with constant coefficients. 

Equation (3.1) can actually be solved explicitly, see, e.g., [18, .30) (Mj- 
A special solution is given by u(t) = e~ Xt , as one can immediately check, 
but the general solution is also available: 


u(t) = e Xt { c\+c 2 


~^e Xv dv 


ci and C 2 being two arbitrary constants. If C 2 / 0,we can see that 

u(t) ~ C 2 T- t ~ v , as t -» +oo, 

A 


for > 0, being A > 0, that is, u(t) decays monotonically as t gets large. 
The oscillatory behavior of the general solution can be seen proceeding 
as follows. Setting u := p(t)v with p(t) = t~ v ^ 2 e~ xt ^ 2 , v will satisfy the 
equation 

v" + q(t)v = 0, (3-2) 


where 


q(t ) : = 


77(2 — 77 ) Xp A 2 

At 2 + 2t ~ T' 


(3.3) 
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Now, as is well known, where q(t) > k 2 > 0, for any arbitrary but fixed 
constant k, the solution is oscillatory in the interval, (t-(k),t .where 

A 77 ± \/2r]\ 2 + 4fc 2 r?(2 — 77 ) 


t± = 

For small k 's we have 

t + (k) — t_{k) 


A 2 + 4 k 2 


A 


k 2 


and hence about t + (k) — t-(k ) ~ 2y / 2r//A. By the Sturm’s comparison 
theorem, the distance between consecutive zeros of u(t ) is less than it/ k. In 
fact, by such a theorem, between any two consecutive zeros of any nontrivial 
solution to u" + k 2 u = 0 , there exists at least one zero of the solution to 
(3.2)-(ph3]). Since i+(fc) —t-(k) (for small k' s) decreases roughly according 


to y^/A, we infer that several oscillations occur in the interval (f_,f_(_). 

On the other hand, the numerical results below show the “physical” 
role of the fractional order of derivatives, a, which is to introduce a rather 
strong oscillatory behavior in the solutions, for some special combinations 
of the parameters A and 77 . 

Generally speaking, varying the parameters A, 77 , and a, in analogy 
with the well-known solutions of the classical damped oscillator equation, 
the solutions can be either overdamped (i.e., they may have a purely mono¬ 
tonic decay), or underdamped (i.e., they may exhibit an oscillatory decay). 
In our numerical solutions, we observe both behaviors, according to the 


physical parameters. In Fig. 4l we plotted the numerical solution of (2.8) 


for several values of A, keeping fixed the other coefficients. We observed an 
underdamped-like oscillatory decay. In Fig. [3j we consider the solutions for 
several values of a and A, for a fixed value of 77 . In this case, for a fixed a, 
changing A, the behavior of the solution switches from an overdamped-like 
monotonic decay to an oscillatory decay pattern, for increasing values of A, 
as expected. A similar less pronounced switch between different behaviors 
can be observed in Fig. [3j varying 77 and a. 

In Fig. [ 6 j we plotted \u'/u\ vs t , for A = 1, 2, 3, correspondingly to the 
numerical solution of ( | 2 . 8 | ). 

In Fig. [t| we plotted \u'/u\ vs t, for A = 5 and large t, correspondingly 
to the numerical solution of ( 2 . 8 ). 
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Figure 4. Numerical solution of (2.8) for a = 1.0,7/ = 
0.5, n = 0.5,/? = 0.5,7 = — 0.5, and different values of A, 
initial discretization step h = 0.1, computed for t E [0,40]. 
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Figure 5. Numerical solution of (2.8) for a = 1.0, 7 = 
0.5,/r = 0.5,/3 = 0 . 5,7 = —0.5, and different values of A, 
initial discretization step h = 0.1, computed for t e [0,10]. 



Figure 6 . Numerical solution of (2.8) plotted as \u'/u\ for 
a = 1.0, 7] = 0.5,// = 0.5, /3 = 0.5,7 = —0.5, and several 
values of A, initial discretization step h = 0.1, computed for 
t G [2,10], 
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Figure 7. Numerical solution of (2.8), plotted as \u'/u\, 
for a = 0.9, A = 5,?? = 0.5, fj, = 0.6, /3 = 0.6,7 = — 0.4, the 
initial discretization step h = 0.1, computed for t E [0,150]. 
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t t 



Figure 8. Numerical solution of (2.8), plotted for several 
values of a E (0,1] and A, with fixed 77 = 0.5, initial dis¬ 
cretization step h = 0.001, computed for t E [0, 5]. 
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a =0.1 


a =0.25 



a. =0.9 



a=1 



Figure 9. Numerical solution of (2.8), plotted for differ¬ 
ent values of a £ (0,1] and 77 , with fixed A = 10, initial 
discretization step h = 0.001, computed for t € [0,5]. 
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